Primeros pasos

El primer paso para realizar el análisis es cargar todas las librerías que utilizaremos y posteriormente cargar la base de datos con las observaciones generadas por el sistema de gestión de flota. Adicionalmente se carga un csv que contiene las coordenadas de todos los paraderos de la ruta Metropolitana.

library(tidyverse)
library(sf)
library(lubridate)
library(purrr)
library(kableExtra)
library(leaflet)
library(leaflet.extras)
library(echarts4r)
observaciones = read_csv("../tracking.csv")
paraderos = read_csv("../02_Metro_paradas.csv")
Tabla de las observaciones del Sistema de Gestión de Flotilla para la ruta Metropolitana
event_date unit_id route_coloquial_name lat lng
19/10/2021 00:19 753 Circuito Metropolitano 21.02143 -89.66084
19/10/2021 00:20 753 Circuito Metropolitano 21.02143 -89.66084
19/10/2021 01:18 753 Circuito Metropolitano 21.02143 -89.66084
19/10/2021 02:17 753 Circuito Metropolitano 21.02143 -89.66084
19/10/2021 02:18 753 Circuito Metropolitano 21.02143 -89.66084
Tabla los paraderos de la ruta Metropolitana
Name title x y
10 10 -89.66343 20.96487
11 11 -89.66343 20.96746
12 12 -89.66344 20.97000
13 SORIANA AV. 128 X AV. MADERO 13 SORIANA AV. 128 X AV. MADERO -89.66345 20.97234
14 14 -89.66302 20.97529

Transformación de clase Sf

Debido a que los data frames anteriores contienen dos columnas para representar la latitud y longitud del punto en cuestión, se puede crear un objeto de tipo sf (simple feature) que permita realizar diversas operaciones espaciales mediante st_as_sf(). Esta agrega una columna adicional al data frame llamada por lo generalmente como geometry que en sí es una lista que contiene las coordenadas del registro. Es necesario declarar el nombre de las columnas en las cuales se encuentran las coordenadas.

observaciones_sf = observaciones %>%
  st_as_sf(coords = c("lng", "lat"))

paraderos_sf = paraderos %>%
  st_as_sf(coords = c("x", "y"))
## Simple feature collection with 609747 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -89.72536 ymin: 0 xmax: 0 ymax: 21.0579
## CRS:           NA
## # A tibble: 609,747 x 4
##    event_date       unit_id route_coloquial_name               geometry
##    <chr>              <dbl> <chr>                               <POINT>
##  1 19/10/2021 00:19     753 Circuito Metropolitano (-89.66084 21.02143)
##  2 19/10/2021 00:20     753 Circuito Metropolitano (-89.66084 21.02143)
##  3 19/10/2021 01:18     753 Circuito Metropolitano (-89.66084 21.02143)
##  4 19/10/2021 02:17     753 Circuito Metropolitano (-89.66084 21.02143)
##  5 19/10/2021 02:18     753 Circuito Metropolitano (-89.66084 21.02143)
##  6 19/10/2021 02:32     753 Circuito Metropolitano (-89.66084 21.02143)
##  7 19/10/2021 02:43     753 Circuito Metropolitano (-89.66084 21.02143)
##  8 19/10/2021 03:16     753 Circuito Metropolitano (-89.66084 21.02143)
##  9 19/10/2021 04:30     753 Circuito Metropolitano (-89.66084 21.02143)
## 10 19/10/2021 05:29     753 Circuito Metropolitano (-89.66084 21.02143)
## # ... with 609,737 more rows

Un problema que emerge del solo transformar la clase del objeto, es que el CRS (Coordinate Reference System) no se encuentra especificado, lo cual puede llevar a operaciones errones. Para asignar un CRS de manera manual al objeto, se recurre a la función st_set_crs() que recibe un objeto de clase sf y un tipo de proyección basado en el sistema EPSG.

observaciones_sf_crs = observaciones_sf %>%
  st_set_crs(value = 4326)

paraderos_sf_crs = paraderos_sf %>%
  st_set_crs(value = 4326)
## Simple feature collection with 609747 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -89.72536 ymin: 0 xmax: 0 ymax: 21.0579
## Geodetic CRS:  WGS 84
## # A tibble: 609,747 x 4
##    event_date       unit_id route_coloquial_name               geometry
##  * <chr>              <dbl> <chr>                           <POINT [°]>
##  1 19/10/2021 00:19     753 Circuito Metropolitano (-89.66084 21.02143)
##  2 19/10/2021 00:20     753 Circuito Metropolitano (-89.66084 21.02143)
##  3 19/10/2021 01:18     753 Circuito Metropolitano (-89.66084 21.02143)
##  4 19/10/2021 02:17     753 Circuito Metropolitano (-89.66084 21.02143)
##  5 19/10/2021 02:18     753 Circuito Metropolitano (-89.66084 21.02143)
##  6 19/10/2021 02:32     753 Circuito Metropolitano (-89.66084 21.02143)
##  7 19/10/2021 02:43     753 Circuito Metropolitano (-89.66084 21.02143)
##  8 19/10/2021 03:16     753 Circuito Metropolitano (-89.66084 21.02143)
##  9 19/10/2021 04:30     753 Circuito Metropolitano (-89.66084 21.02143)
## 10 19/10/2021 05:29     753 Circuito Metropolitano (-89.66084 21.02143)
## # ... with 609,737 more rows

Ahora vemos que nuestros objetos tienen asociado el CRS WGS 84.

Creación de buffers

El siguiente paso para nuestro análisis, es crear un buffer para los puntos de los paraderos. La distancia de estos buffers deberá ser de 20 mts. Para esto, se utilizará la función st_buffer() que recibe un objeto de tipo sf y la distancia del buffer.

Sin embargo, ya que que el CRS WGS 84 se basa en grados, es necesario antes realizar una transformación a otro CRS diseñado con base en metros. En ese sentido, se utilizará el EPSG: 6372 Mexico ITRF2008 / LCC utilizando la función st_transform().

paraderos_buffer = paraderos_sf_crs %>% 
  st_transform(crs = 6372) %>%
  select(Name, geometry) %>%
  st_buffer(dist = 20)

Intersección de observaciones con buffer

Debido a limitaciones del sistema, no es posible conocer el momento exacto en el que una unidad pasó por un paradero. Para solucionar esto, se definió un buffer de 20 mts en cada paradero y se filtrarán únicamente aquellas observaciones que se encuentren dentro de alguno de estos. Para lograr esto se utiliza la función st_filter() que recibe dos objetos de tipo sf y aplica una función de geométrica, en este caso, la intersección st_intersects.

Se debe mencionar que como la proyección del buffer usa el CRS EPSG: 6372, se realizó una reproyección para las observaciones.

observaciones_repro = observaciones_sf_crs %>%
  st_transform(crs = 6372)

observaciones_filtered = st_filter(observaciones_repro, paraderos_buffer, .pred = st_intersects)

Asociación de observaciones con paraderos

Debido a que existen paraderos cuyos buffers se traslapan entre sí, en muchas ocasiones, una misma observación puede estar contenida dentro de dos o más paraderos, por tanto, se asociará cada observación a el paradero más cercano, evitando así duplicidades en los registros. Lo anterior se realiza con la función st_join() que asocia un paradero a cada observación; la función st_nearest_feature() modifica el comportamiento de la primra función, asegurándose que a cada observación se le asocie al paradero más cercano.

observaciones_nearest = st_join(observaciones_filtered, paraderos_buffer, join = st_nearest_feature)
Observaciones depuradas
event_date unit_id route_coloquial_name geometry Name
19/10/2021 10:13 753 Circuito Metropolitano POINT (3782833 1047811) 275 UNIDAD DEPORTIVA KUKULKAN-CBTIS 95
19/10/2021 10:53 753 Circuito Metropolitano POINT (3784925 1053110) 328
19/10/2021 10:57 753 Circuito Metropolitano POINT (3784674 1054691) 335
19/10/2021 12:15 753 Circuito Metropolitano POINT (3776638 1053267) 433 PLAZA DORADA
19/10/2021 13:14 753 Circuito Metropolitano POINT (3782804 1047744) 183 CBTIS 95
19/10/2021 16:26 753 Circuito Metropolitano POINT (3784112 1048831) 172 UNIVERSIDAD PEDAGÓGICA
19/10/2021 16:29 753 Circuito Metropolitano POINT (3784042 1050012) 291
19/10/2021 17:24 753 Circuito Metropolitano POINT (3781877 1056628) 362 BOXITO KALIA
19/10/2021 18:02 753 Circuito Metropolitano POINT (3776760 1057372) 411 GLORIETA LA MESTIZA
19/10/2021 18:24 753 Circuito Metropolitano POINT (3776225 1053368) 432 WALMART PENSIONES

Procesamiento de datos

Una vez que ya se cuenta con una base de datos pre-procesada de las observaciones que se encuentran lo suficientemente cerca de un paradero, el paso siguiente es aplicar un procesamiento más detallado para la final presentación de los hallazgos.

En este punto, los datos sobre las coordenadas de las observaciones ya no son imperativas, de forma que son removidas del análisi para hacer las operaciones más ágiles.

observaciones_pre = observaciones_nearest %>%
  tibble() %>%
  mutate(
    event_date = dmy_hm(event_date)
  ) %>%
  select(-geometry, -route_coloquial_name)

Limpieza de datos

Debido a que por cada minuto las unidades envían información acerca de su ubicación, el número de pasajeros, etc., es muy probable que se cuenten con observaciones de una misma unidad que se mantuvo durante más de un minuto en un mismo paradero. Estas observaciones tendrán que ser filtradas del análisis.

Los datos primero se ordenan por el nombre del paradero, luego por la fecha del posteo. Posteriormente, se agrupan de acuerdo a la fecha del posteo, del nombre del paradero y del ID de la unidad. Esto permite observar los momentos en los que la unidad x pasó por el paradero y en el día z.

observaciones_nested = observaciones_pre %>%
  arrange(Name, event_date) %>%
  group_by(date = as_date(event_date), Name, unit_id) %>%
  nest()
## # A tibble: 57,409 x 4
## # Groups:   unit_id, Name, date [57,409]
##    unit_id Name           date       data            
##      <dbl> <chr>          <date>     <list>          
##  1     800 1 talleres ado 2021-10-08 <tibble [4 x 1]>
##  2     711 1 talleres ado 2021-10-08 <tibble [1 x 1]>
##  3     743 1 talleres ado 2021-10-08 <tibble [2 x 1]>
##  4    1093 1 talleres ado 2021-10-08 <tibble [1 x 1]>
##  5     731 1 talleres ado 2021-10-08 <tibble [1 x 1]>
##  6     706 1 talleres ado 2021-10-08 <tibble [2 x 1]>
##  7     716 1 talleres ado 2021-10-08 <tibble [1 x 1]>
##  8     726 1 talleres ado 2021-10-08 <tibble [2 x 1]>
##  9     740 1 talleres ado 2021-10-08 <tibble [1 x 1]>
## 10     732 1 talleres ado 2021-10-08 <tibble [1 x 1]>
## # ... with 57,399 more rows

El anterior data frame es un ejemplo de un nesteo, en el cual, una columna de la tabla es una lista de tablas. Con este data frame, se aplicará una función que permite calcular el tiempo que le tomó a la unidad volver a pasar por cada paradero, lo cual finalmente nos permitirá filtrar aquellas observaciones en donde este tiempo sea menor a 5 minutos.

Esta función es una función propia con nombre timediff_na() la cual con base en una columna de fechas, ordenada previamente de más antigua a más reciente calcula el tiempo en minutos entre observaciones consecutivas.

timediff_na = function(df) {
  minutos = interval(start = lag(df$event_date, n = 1), end = df$event_date,) %/% minutes(1)
  
  return(minutos)
}

Esta función se mapea mediante la función map() de la librería purrr la cual aplica una función determinada a una lista y regresa también una lista.

observaciones_map = observaciones_nested %>%
  mutate(
    minutos = map(data, timediff_na)
  )

Posteriormente, se filtran aquellos valores que están en un rango de 1 a 5 minutos, debido a que esto querría decir que a la unidad le tomó entre 1 a 5 minutos dar una vuelta al circuito, lo cual no es factible.

observaciones_filtrado = observaciones_map %>%
  unnest(c(data, minutos)) %>%
  replace_na(list(minutos = 10)) %>%
  filter(!between(minutos, 1, 5))

Cálculo de minutos entre unidades por paradero

Una vez depurado el data frame, se procede a realizar un procedimiento análogo, primero nesteando el data frame y luego calculando los minutos entre arribos.

observaciones_final = observaciones_filtrado %>%
  ungroup() %>%
  select(-minutos, -unit_id) %>%
  arrange(Name, event_date) %>%
  group_by(date, Name) %>%
  nest() %>%
  mutate(
    minutos = map(data, timediff_na)
  ) %>%
  unnest(c(data, minutos)) %>%
  select(event_date, date, Name, minutos)
Minutos entre arribo por paradero
event_date date Name minutos
2021-10-08 06:12:00 2021-10-08 1 talleres ado NA
2021-10-08 06:29:00 2021-10-08 1 talleres ado 17
2021-10-08 07:31:00 2021-10-08 1 talleres ado 62
2021-10-08 07:55:00 2021-10-08 1 talleres ado 24
2021-10-08 08:14:00 2021-10-08 1 talleres ado 19

La observación del primer arribo del día por paradero siempre quedará indeterminado, puesto que no hay ningún arribo con el cual compararlo.

Presentación de hallazgos

Con la información ya procesada y debidamente depurada, se estima que para la ruta Metropolitana, durante el 08 de octubre al 22 de octubre, el tiempo de espera promedio fue de 64 minutos.